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ABSTRACT 

Parametrisation of galaxy morphologies is a challenging task, for instances in shear 
measurements of weak gravitational lensing or investigations of formation and evolu- 
tion of galaxies. The huge variety of different morphologies requires a parametrisation 
scheme that is highly flexible and that accounts for certain morphological observables, 
such as ellipticity, steepness of the radial light profile and azimuthal structure. In this 
article, we revisit the method of sersiclets, where galaxy morphologies are decomposed 
into a set of polar basis functions that are based on the Sersic profile. This approach 
is justified by the fact that the Sersic profile is the first-order Taylor expansion of any 
real light profile. We show that sersiclets indeed overcome the modelling failures of 
shapelets in the case of early-type galaxies. However, sersiclets implicate an unphysical 
relation between the steepness of the light profile and the spatial scale of the polyno- 
mial oscillations, which is not necessarily obeyed by real galaxy morphologies and can 
therefore give rise to modelling failures. Moreover, we demonstrate that sersiclets are 
prone to undersampling, which restricts sersiclet modelling to highly resolved galaxy 
images. Analysing data from the weak-lensing GreatOS challenge, we demonstrate 
that sersiclets should not be used in weak-lensing studies. We conclude that although 
the sersiclet approach appears very promising at first glance, it suffers from conceptual 
and practical problems that severly limit its usefulness. In particular, sersiclets do not 
provide high precision results in weak-lensing studies. Finally, we show that the Sersic 
profile can be enhanced by higher-order terms in the Taylor expansion, which can 
drastically improve model reconstructions of galaxy images. When orthonormalised, 
these higher-order profiles can overcome the problems of sersiclets while preserving 
their mathematical justification. However, this method is computationally expensive. 

Key words: Galaxies: general - Methods: data analysis, statistical - Techniques: 
image processing ~ Gravitational lensing: weak. 



1 INTRODUCTION 

There are two scientific key diagnostics in modern inves- 
tigations of galaxy formation and evolution, namely star- 
formation rates and galaxy morphologies. Galaxy morpholo- 
gies are known to correlate to different degrees with other 
more physical parameters such as star formation (e.g. 



Ken- 



nicutt|1998j ) or environment (e.g., van der Wei et al.|2010 | 



Although morphology by itself is not a fundamental pa- 
rameter of a galaxy, it provides a direct observable. Con- 
versely, e.g., star-formation rates are derived quantities 
based on additional assumptions and extrapolations (e.g., 
Rosa-Gonzalez et al. 2002 J . Studies of galaxy morpholo- 



gies are therefore an important complementary means for 
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studying the physics of galaxy formation. Furthermore, cer- 
tain investigations are solely based on morphologies, e.g., 
weak-lensing measurements (e.g., Bernstein fc Jarvis||2002 1 
or studies of angular-momentum alignments of spiral galax- 



ies (e.g., Slosar et al. 2009 Andrae & Jahnke 



m prep. 



Concerning weak gravitational lensing, the investigation of 
shape-measurement algorithms currently is a very active and 
important field of research since established and well known 
procedures fail to measure object shapes accurately enough 
to fully exploit the potential of future large-scale imaging 



surveys (e.g. Bridle et al. 2010 Bernstein 20101. Concern- 



ing morphological classification, galaxies were traditionally 
described mostly by visual classifications. However, given 
the enormous amount of galaxies known from today's sky 
surveys, it is evident that an automated parametrisation 
of galaxy morphologies is required. Consequently, there has 
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been substantial eflort to define automated parametrisation 



schemes for galaxy morphologies (e.g., 


Sersic 1968 


Abraham 


et al. l996, 2003, Bcrshady et al.|2000 


|Conselice'2003||Lotz 



et al. 2004, 2008, to name just a few). Unfortunately, these 



parametrisation schemes usually invoke rather restrictive as- 
sumptions, such that they lack the flexibility to describe 
the huge variety of different galaxy morphologies present in 
modern databases (Andrae et al. 2010a I. The Sersic pro- 



file ( jSersic 1968i ) is one example in this respect: The Sersic 
model is the first-order Taylor expansion of any real radial 



light profile (cf. Appendix Al I and therefore provides a good 
first-order approximation to real galaxies. However, second- 
order effects - e.g., deviations of the radial profile from the 
Sersic form or the appearance of azimuthal structures such 
as star-forming regions, galactic bars or spiral-arm patterns 
- are typically not negligible, and hence there have been 
several attempts to enhance the Sersic profile: 



• Mixtures of two or more Sersic profiles (e.g. [Simmat] 
et al.|2010| ). 

• Modification of the Sersic profile by Fourier modes in 



order to describe azimuthal structures (e.g., Galfit 3, Peng 
et al.|2010| |. 



Orthogonalisation of the Sersic profile in order to de- 
scribe azimuthal structures (Ngan et al.||2009l. 
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Figure 1. Radial basis functions if orthonormalised on the inte- 
gration interval r £ [0, 0\ using a Sersic profile with index ng = 2 
as weight function | |Ngan e t al."2009|. The vertical dashed fines 
indicate r = fi, where the orthonormalisation interval ends. The 
Sersic profile is not able to suppress the polynomials outside the 
orthonormalisation interval. 



2 POLAR SERSICLETS 

We briefly comment on a prior attempt to introduce 
sersiclets before we then give our definition. 



The idea of using basis functions in order to parametrise 
galaxy morphologies is not new. For instances, "shapelets" 
pjfrcgicr 2003| are a set of basis functions based on the 
Gaussian profile. The appealing properties of shapelets gave 
rise to a number of investigations of galaxy morphologies 



(e.g., Kelly & McKay|2004 


2005 Andrae et al. 


2010b 1 and 


weak-lensing studies (e.g., 


Massey et al. 2007 


Ferry et al. 


2008). However, as was recently shown by Melchior et al. 



1 2010 1, shapelets can suffer from strong biases that originate 
from the radial light profiles of galaxies potentially being 
much steeper than a Gaussian profile, in particular for early- 
type galaxies. These biases manifest themselves as ring-like 



artefacts in the models and the residual maps (cf. Sect. 3.2 1. 
They are able to wash out virtually all the information about 
ellipticity of an object and are likely to affect other morpho- 
logical quantities in a similar way. Hence, it was obvious to 
introduce a set of basis functions based on the Sersic profile 
( Ngan et al.||2009 l, called "sersiclets". As these basis func- 
tions explicitly account for the steepness of radial profiles, 
they are highly fiexible and seem to be promising candi- 
dates to provide a method for describing the huge variety of 
different galaxy morphologies. 

This article is organised as follows. We define sersiclet 
basis functions in Sect. [2] Section[3]presents systematic tests 
of sersiclets, working out their reliability and limitations. We 
investigate the performance of sersiclets in an application to 
artificial weak-lensing data in Sect. |4] In Sect. [5] we discuss 
the potential of enhancing the Sersic profile by higher-order 
Taylor expansions and briefly investigate the orthonormali- 
sations of such light proflles. We conclude in Sect. [6] Further 
details are given in the appendices. Appendixj^contains the 
proof that Sersic profiles are first-order Taylor expansions. 
Appendix |B] provides the analytic derivation of the sersiclet 
basis functions. Finally, we explain in Appendix [C] how to 
fit sersiclet models to imaging data. 



2.1 The first attempt to introduce sersiclets 

Given the aforementioned problems of shapelets, it was ob- 
vious to orthonormalise the Sersic profile. Moreover, this 
choice is "natural" because the Sersic profile is the first-order 



Taylor expansion of any real light profile (cf. Appendix Al I. 
This is likely to remove the strong biases of shapelets. The 
resulting basis functions are called sersiclets and were first 



investigated by Ngan et al. ( 2(]09| . 

However, Ngan et al. ( 2009 1 faced a severe problem with 



sersiclets. In a simple test case, Ngan et al. ( 2009 1 observed 
that the polar sersiclets were incapable of fitting a given 
object. Even when increasing the number of basis functions 
used for the decomposition, the model did not converge to 
the given image data. This directly implies that the basis 



functions constructed by Ngan et al. ( 2009 1 are not com- 



plete. We found out that this problem originates from the 



fact that Ngan et al. ( 2009 1 orthonormalised their basis 
functions out to one half-light radius - not on the interval 
r G [0, oo[ - but then extended their basis functions beyond 
this range. Figure [T] shows what happens in this case. In the 
inner region, the basis functions may be orthonormal, but for 
larger radii the leading order monomial in the basis functions 
dominates over the Sersic profile. This creates an artificial 
bump outside the region of orthonormality, which appears 
in all radial modes of order I > and becomes the dominant 
feature. Consequently, these basis functions loose their linear 
independence, which explains the non-convergence observed 



by Ngan et al. (20091 



In order to solve this problem, they suggested to discard 
all modes with azimuthal "quantum number" m 7^ 0, i.e., all 
mode with azimuthal structure, and only maintain the radial 
modes. While this approach may be viable for weak-lensing 
applications, where azimuthal structures are rarely visible, 
it does not solve the actual problem, which is the extension 
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Figure 2. Cartesian ground state using Gaussian profile (left 
panel) and Sersic profile with index ng = 2 (right panel). Unless 
we employ the Gaussian profile, the Cartesian weight function is 
not spherically symmetric. 



beyond the orthonormalisation interval. Moreover, the re- 
sulting basis functions may still represent an expansion into 
radial modes, but it is not only radial but also azimuthal 
structure of galaxies that is interesting, e.g., in morpho- 
logical classification. In fact, it is the ability to describe 
azimuthal structure in a (theoretically) well defined way 
that makes basis-function expansions of galaxy morpholo- 
gies such a compelling approach. Consequently, [Jimenez^ 
Teja & Bemtez (20111 criticised that this implementation 



of sersiclets is incapable of describing irregular galaxies. 



2.2 Definition of polar sersiclets 

We now define the sersiclet basis functions. In the case of 
shapelets, we can define both Cartesian and polar shapelets 
due to the mathematical features of the Caussian weight 
function. However, for general Sersic profiles, we cannot 
meaningfully define Cartesian basis functions because the 
ground states hardly resemble any known galaxy morphol- 
ogy, as is shown in Fig. [2] Therefore, we have to introduce 
sersiclets as polar basis functions. 

The crucial idea to overcome the problems faced by 
Ngan et al. ( 2009 1 is to realise that we cannot expect a 



galaxy's structure being well captured within one half-light 
radius. For instances, disc galaxies may exhibit spiral-arm 
patterns in their outskirts. Therefore, we require orthonor- 
mality on the full interval r € [0, oo[. Apart from solving the 
problems reported by Ngan et al.| ( 2009 1, this also has the 



advantage that the orthonormal polynomials exist analyti- 
cally. The sersiclet basis functions are derived in Appendix [B] 
and are given by 



cxp 



1/ns 



(1) 



where the normalisation factor is also derived in Ap- 
pendix [B] Apart from the linear expansion coefficients, 
sersiclets have the following nonlinear model parameters: 

• The maximum order A'^max- 

• The Sersic index ns- 

• The scale radius /3. 

• The centroid position xq, which enters via r = \x — xq\. 

• The complex ellipticity e, which enters via a shear trans- 



formation (e.g. Bartelmann & Schneider 2001 1 thereby af- 
fecting 



^ It is advantageous to use the real and imaginary parts of the 
complex ellipticity as model parameters, rather than orientation 



The coefficient b is defined by ( Graham fc Driver|2005 1 

r(2ns) = 27(2ns,fe), (2) 

where F and 7 denote the complete and incomplete gamma 
functions. We explain in Appendix[C]how to fit these model 
parameters for given imaging data. In the case of ns ~ 0.5 
and b = 1 the polar sersiclets reduce to the special case of 
polar shapelets]^ Furthermore, for ns ~ 1 and 6=1, we 
get a set of basis function that could be called "disclets", 
since they have an exponential profile as weight function. 
Similarly, we could define "deVaucouleurlets" . All these sets 
are special cases of Eq. ([TJ . Finally, we emphasise that when 
fitting a sersiclet model, we always also fit a Sersic model, 
which is the ground state of the sersiclets. Figure |3] displays 
an example of sersiclet basis functions. 



2.3 Interpretation of the Sersic index 

In the case of sersiclet basis functions, the Sersic index ns 
changes its interpretation. First, ns regulates the steepness 
of the weight function, which is a normal Sersic profile. Sec- 
ond, via the steepness of the weight function, ns also reg- 
ulates the spatial scale on which the associated Laguerre 
polynomials oscillate. In simple words, there is a fixed rela- 
tion between steepness of the weight function and oscillation 
scale of the polynomials. However, real morphologies do not 
necessarily obey such a relation. For instances, the steepness 
of the bulge is not related to the positions of spiral arms or 
star-forming regions. In practice, this can lead to modelling 
problems, if the steepness of the profile and the range of os- 
cillation scales required for a faithful description of a given 
galaxy morphology do not match. 



3 TESTING SERSICLETS 

Before we can apply sersiclets to scientific questions, we need 
to investigate the performance and fidelity of sersiclets. 



3.1 Completeness 

First, we need to test the completeness of sersiclets in or- 
der to proof that we indeed overcome the problems reported 



by Ngan et al. (20091. We decompose three real galaxy im- 



ages from the Sloan Digital Sky Survey into sersiclets with 
increasing maximum order. If our basis functions were not 
linearly independent, completeness would break down and 
the x^-values would not decrease with increasing maximum 
order as in the case of |Ngan et al.| (|2009|. Figure [4] shows 
the results. First and foremost, the residuals of sersiclets are 
decreasing with increasing maximum order, i.e., we indeed 
set up a set of basis functions that are linearly indepen- 
dent. Furthermore, it is interesting to compare the resid- 
uals of sersiclets with those of circular shapelets. In the 



angle and axis ratio. Using the orientation angle as parameter 
would cause severe problems with convergence in the case of 

nearly spherical objects. 

^ Note that the polar functions defined by [Massey 



Refregier] 



l j2005| are not orthogonal. A correct functional form of polar 
shapelets is derived in Appendix [B] 
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Figure 3. Polar scrsiclct basis functions with ng = 1, 6 = 1, and 
e = 0. The real components of the complex-valued basis functions 
are shown in the top panel, the imaginary components in the 
bottom panel. The basis functions with m = are wholly real. We 
note that the polar basis functions exhibit lots of substructure in 
the central region. Moreover, these central substructures become 
very small with increasing radial order I. 



case of the spiral galaxy, sersiclets yield residuals compa- 
rable to shapelets. This is not surprising, because shapelets 
excel in modelling extended objects with lots of substruc- 
ture, such as a face-on spiral galaxy. However, in the case of 
the edge-on disc and the elliptical galaxy, the x'^ -values of 
sersiclets are substantially lower. Evidently, sersiclets out- 
perform shapelets in modelling galaxies of these types, as 
expected from the steep light profiles. 
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Figure 4. Dependence of of sersiclet decomposition (solid 
lines) and shapelet decompositions (dashed lines) on maximum 
order, Afmax, for the spiral galaxy, the edge-on disc, and the el- 
liptical galaxy shown in Fig. [s] For each example, we also show 
the number of pixels as a horizontal dotted line. 



3.2 Image decompositions 

We have seen in the previous section that sersiclet decom- 
positions produce substantially lower residuals than (circu- 
lar) shapelet decompositions when it comes to modelling 
galaxies that exhibit steep profiles, such as elliptical galax- 
ies or edge-on discs. However, we still have to check whether 
sersiclets indeed overcome the ring-like artefacts produced 
by shapelets. 

Figure [5] compares the best-fitting models and resid- 
uals of (circular) shapelets and (elliptical) sersiclets using 
A^'max = 8 for the three test galaxies. As expected from 
the similar x^-values, in the case of the spiral galaxy, both 
shapelets and sersiclets perform well in modelling the spiral- 
arm patterns. Sersiclets fit the central region better than 
shapelets, whereas shapelets tend to describe the outskirts 
better. As expected from a sersiclet model with ns > 0.5, 
the polynomial oscillations appear on smaller scales than for 
the shapelet model with ng = 0.5. In the case of the edge-on 
disc, sersiclets are clearly superior. First, the ring-like arte- 
facts of shapelets are gone. Second, the intrinsic ellipticity of 
the sersiclet model allows the basis functions to also describe 
the outermost regions of the disc, whereas these regions go 
almost unfitted by circular shapelets. In the case of the el- 
liptical galaxy, the ring-like artefacts are very prominent in 
the shapelet residuals. Conversely, the sersiclet residuals do 
not exhibit any artefacts of this kind, i.e., the steepness of 
the Ught profile is indeed described properly. 



3.3 Orthonormality and sampling 

Let E denote the pixel covariance matrix of the noise in the 
imaging data and X the design matrix (see Appendix |C1[ ). 
In order to get the maximum-likelihood estimate of the ex- 
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shapelet residuals shapelet model data sersiclet model sersiclet residuals 




Figure 5. Models and residual maps resulting from sersiclets and shapelets for the spiral galaxy (top row), the edge-on disc (central 
row), and the elliptical galaxy (bottom row). All models used A'^max = 8. For the sake of visualisation, the colour code of the model and 
data maps is nonlinear and ranges from — 4.5(T to the maximum value. The colour code of the residual maps is linear and ranges from 
— 4.5iT to +4.50", where cr denotes the standard deviation of the background noise. The peak significance of the imaging data is 132(T for 
the spiral galaxy, lOlo" for the edge-on disc, and 610o- for the elliptical galaxy. 



pansion coefficients from Eq. (C3l, the matrix X"^ • • X 
- or ■ X for uncorrelated pixel noise - needs to be in- 
vertible. If it were not for pixellation, finite image grids, and 
the PSF, X^ ■ X should be an identity matrix due to the 
orthonormality of the basis functions. We now investigate 
the orthonormality of sersiclets in two tests. 

First, we compare the orthonormality of polar shapelets 
(sersiclets with ns = 0.5 and 6=1) and Cartesian shapelets 
of [Melchior et al. (20071. We sample polar and Cartesian 
shapelets of constant scale radius /3 — 5 on pixel grids of 
sizes 50x50, 100x100, and 500x500 and compute the de- 
terminant of X^ ■ X for maximum orders A'max ^ 16, 
where det(X'^ • X) < 1 indicates a violation of orthonormal- 
ity. Figure [6] shows the results of this test. In the case of the 
50x50 grid (panel (a)) both Cartesian and polar shapelets 
suffer from non-orthonormality for increasing A^maxQThis 
problem can be caused either by boundary truncation or 
undersampling of higher-order modes (with rapidly oscillat- 



^ In particular, for polar shapelets and A^max ^ 8 the determinant 
of ■ X is zero, i.e., X^ ■ X is not invertible and the estimator 



of the expansion coefficients (Eq. (C3i) breaks down 



ing polynomials). Therefore, in panel (b) we increase the 
pixel grid from 50x50 to 100x100. The non-orthonormality 
of Cartesian shapelets is now cured, i.e., it has indeed been 
caused by boundary truncation. However, polar shapelets 
still exhibit non-orthonormality. Therefore, in panel (c) we 
increase the pixel grid again now from 100 x 100 to 500x500. 
The behaviour of polar shapelets is unchanged, i.e., the non- 
orthonormality is not caused by boundary trunction. In or- 
der to demonstrate that the non-orthonormality is caused 
by undersampling, we increase in panel (d) the object size 
from /3 = 5 to /3 = 25 while keeping the 500x500 pixel 
grid. The non-orthonormality of polar shapelets is now al- 
most cured, which verifies that this was an undersampling 
effect. Evidently, undersampling has a larger impact on po- 
lar shapelets than on Cartesian shapelets. The reason is that 
polar basis functions exhibit most of their structure in the 
central region (cf. Fig. [sj and therefore suffer strongly from 
pixellation. Loosely speaking, polar basis functions do not 
appreciate being sampled on a Cartesian pixel grid. 

Second, having attested problems of polar shapelets 
with orthonormality due to undersampling, we now investi- 
gate the orthornomality of (polar) sersiclets in general. We 
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Figure 7. Smallest (solid line) and largest (dashed line) diagonal elements of matrix • X as a function of j3 at Mnax = 12. Deviation 
from unity (dotted line) indicates the loss of orthonormality. The image size was 100x100 pixels (left columns) and 1000x1000 pixels 
(right column). We used b = 2 log 3, such that at one scale radius the weight function drops to one third of its central value. 



compare the value of the largest and smallest diagon al ele- 



ment of -X, respectively, repeating the same test as Berry 



|et al.| ( |2004[ ) did for shapelets. We evaluate sersiclet models 
of maximum order A^'max = 12, Sersic indices ns = 0.5, 1, 2, 4 
and axis ratios q = 1,0.8,0.6,0.4 on a 100x100 pixel grid 
for varying scale radii p. There is no PSF in this test. Fig- 
ure [7] shows test results. For the moment, we only consider 
the first row in Fig. [7| which corresponds to polar shapelets. 
The results agree with Fig. [6] revealing undersampling ef- 
fects for small /3 and boundary truncation for large /3 (cf. 
Melchior et al.|[2007 and discussion therein). Furthermore, 
we can see that the axis ratio has only a mild impact on the 
orthonormality. However, inspecting the other rows corre- 
sponding to steeper profiles (ns = 1,2,4), Fig. [t] reveals se- 
rious violations of orthonormality. Obviously, for ns ^ 1 the 
undersampling regime is not overcome before the boundary 
effects set in. This is confirmed by the right-most column in 
Fig. [T] which used an enlarged pixel grid while keeping the 
resolution (scale radius /3) constant. In this case, sersiclets 
with ns = 1 exhibit decent orthonormality on this larger 
grid, while the cases with ng > 1 still to not reach an ac- 
ceptable level of orthonormality. The reason for this pecu- 
liar behaviour is the argument h (r/P)^^"^ of the associated 



Laguerre polynomials in Eq. ([T|. It implies that the polyno- 
mials are only slowly varying with r for large values of ns. 
Consequently, polar sersiclets have a serious problem with 
orthonormality, especially for large Sersic indices. 

What is the impact of undersampling? If det(X"^ • X) 
is too close to zero, numerical inaccuracies will dominate 
its inversion and the estimate of the expansion coefficients 



(Eq. (C3l) will catch up random errors. Consequently, this 
increases the uncertainty in the model parameters and may 
even lead to completely random parameter values. 

The impact of undersampling could be alleviated by 
sampling the model several times within each pixel. In the 
limit of infinitely fine sampling this amounts to a convolu- 
tion of the model with the pixel-response function. As the 
undersampling problem stems from a key feature of the em- 
ployed sersiclet model, namely rapid polynomial oscillations, 
an oversampled model may better approximate the image 
data. However, the information about the galaxy morphol- 
ogy contained in the expansion coefficients has been lost and 
is not restored by oversampling. Instead, excessive oversam- 
pling generates artificial information that compromises the 
model parameters. As we require meaningful information 
from the expansion coefficients, we will not further pursue 
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Figure 6. Determinant of X"^ ■ X for Cartesian (dashed lines 
with squares) and polar (solid lines with dots) shapelets for in- 
creasing maximum order Afmaxi different scale radii /? and grid 
sizes. Panel a: At large A^max both Cartesian and polar shapelets 
deviate from det(X-^ ■ X) = 1. Panel b: Cartesian shapelets main- 
tain det(A:^ ■ X) = 1 for all Nm ax) proving that boundary ef- 
fects are now negligible. Panel c: Polar shapelets still differ from 
det(X^ ■ X) = 1 at large Amax on this enlarged grid, i.e., this is 
not a boundary effect. Panel d: Polar shapelets now also maintain 
det(X^ ■ X) = 1, i.e., the remaining effect was due to undersam- 
pling. 



this approach. If one wants to employ this approach, one has 
to bear in mind, that the additional convolution will lead to 
covariances among the coefficients and therefore render the 
fitting procedure more complicated. 

3.4 Analogy to Nyquist frequency 

In order to avoid the undersampling problem in practice, 
e.g., in morphological classification, we want to derive a 
lower limit to the scale radius /3 of a sersiclet model. This 
limit is given by comparing the pixel size to the scale on 
which the radial components of sersiclets - essentially the 
Laguerre polynomials - vary. In fact, this is an analogy to 
the Nyquist frequency in the case of Fourier transform. 

The key to set this lower limit is to identify the short- 
est scale on which a Laguerre polynomial varies. This scale 
can be inferred from the roots of the Laguerre polynomial. 
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5.4373 


21.632 


67.371 


180.49 


12 


2.4002 


5.8922 


25.137 


83.305 


236.08 


13 


2.4944 


6.3471 


28.902 


101.54 


303.46 


14 


2.5853 


6.8018 


32.926 


122.23 


384.13 


15 


2.6730 


7.2565 


37.210 


145.53 


479.75 


16 


2.7579 


7.7110 


41.754 


171.59 


592.02 


17 


2.8403 


8.1656 


46.557 


200.57 


722.76 


18 


2.9204 


8.6201 


51.620 


232.61 


873.87 


19 


2.9983 


9.0745 


56.943 


267.88 


1047.4 


20 


3.0742 


9.5289 


62.525 


306.52 


1245.3 



Table 1. Lower limits of scale radii /3 in units of pixels for different 
maximum radial orders Amax according to Eq. Q . Here we used 
b = 2ns - 1/3. 



An associated Laguerre polynomial Li{x) of order I has / 
real-valued, positive roots within the interval {0,1 + k + {I 

+ k]. The smallest scale that can be resolved by this 
polynomial is given by the distance between a; = (the 
peak of the radial component) and the first root xi > 0. 
Unfortunately, the roots of associated Laguerre polynomials 
are not known analytically, so the first root xi has to be 
inferred numerically. As the argument of the associated La- 
guerre polynomial is x = b[r / jSY^"^ , given x\, we can infer 
the radius of the first root to be r\ = I3(xi/b)^^ . Hence, the 
distance between r = and the first root is given by 



Ar- = ri - = ri = /?(a;i/&)"^' 



(3) 



This scale Ar should be well resolved by the pixel grid in 
order to avoid undersampling. An optimistic lower limit is 
Ar = 1, i.e., the radial component drops from its central 
peak to its first root within a single pixel. Setting Ar = 1 
and solving for the minimal /3, we obtain. 



P ^ /3min = {b/xtY 



(4) 



where /3 and hence /3min are given in units of pixels. Table [T] 
provides values of /3min for some realistic values of ns and 
maximum polynomial order A^max = I- Considering this ta- 
ble, we also have to keep in mind that the image grid has to 
be large enough such that several scale radii fit into it. Evi- 
dently, as ns increases, the limit becomes larger and larger. 
This can also be inferred qualitatively from Fig. [7] but Eq. 
(j4]| provides a quantitative result. Obviously, for ns ^ 2 
the lower limit becomes extraordinarily large. In particular, 
71S = 4 would require an extremely well resolved object with 
radius of several thousand pixels. 

In the case of shapelets, it is standard practice to choose 
a maximum order A'^max for a decomposition according to 
signal-to-noise ratio and resolution. For sersiclets, the res- 
olution constraint plays a much more important role and 
restricts models with large to low Nmmx, with which com- 
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Set 


91 


92 


0007 


0.0005405 


0.0069236 


0026 


-0.0527875 


-0.0090224 


0035 


0.0166067 


-0.0045223 


0048 


0.0708862 


-0.0377040 


0056 


0.0325078 


0.0978346 


0091 


-0.0246346 


-0.0488837 


0126 


0.0170977 


-0.1383142 


0135 


0.0596913 


0.0416342 


0268 


-0.0653126 


-0.0883511 


0281 


-0.0431769 


0.0462176 



Table 2. Data sets chosen from the GreatOS sample and their 
complex shears 9 = 9i + ig2- 



plicated morphological features could receive improper de- 
scription. As elliptical galaxies rarely exhibit such features 
though, this restriction may not be overly problematic. 



4 APPLICATION TO WEAK-LENSING DATA 

Given the fact that shapelets were also intended for shear 



measurements in weak-gravitational lensing (e.g., Refregier 
&c Bacon 2003 Chang et al. 20041, we now study the po- 



tential of sersiclets in this field of research. From the sim- 
ulated GreatOS data ( Bridle et al. [2009 ) we selected ten 
images with artificial galaxies, whose shear values are given 
in Table 2j^ Every set contains 10,000 objects sampled on 
a 40x40 pixel grid. Within such a set the applied shear is 
constant and the aim is to retrieve its value. Furthermore, 
all objects have been convolved with a PSF that is a Moffat 
profile with FWHM= 2.85, (3 = 3.5, and intrinsic ellipticity 
e = 0.019 - 0.007i (cf. [Bridle et al.|2010| ). 

We decompose all objects into sersiclets with maximum 
orders A^'max = 0, 2, 4, 6, taking into account the PSF by for- 
ward modelling. For each set, the shear is estimated via the 
mean ellipticity e from the 10,000 artificial galaxies via the 
ellipticity parameter of the sersiclet model, which theoret- 
ically provides an unbiased estimator of the gravitational 
shear (cf. Bartelmann fc Schneider, |2001[ |. Figure [s] shows 
the test results as input vs. estimated shear. The goodness of 
the shear estimates is parametrised by a straight-line model 



f{x) = a + bx . 



(5) 



A perfect shear estimator yields a = and b — 1. If a real 
shear estimator yields an offset a 7^ 0, i.e., a shear is detected 
although the input shear was zero, this typically implies that 
the PSF is not properly corrected for. If 6 < 1, the true 
shears are underestimated. 

Given the test results shown in Fig.jsj the sersiclets per- 
form best for A^ax = 0, where the slope b is consistent with 
1 (at a fairly large offset a). This is not surprising, because 
for Amax = sersiclets reduce to pure Sersic profiles, which 
were used to simulate this subset of the GreatOS datalf] 



This restriction on only ten sets was necessary, because 
analysing the whole GreatOS data set would have been too time- 
consuming. 

^ The data sets were taken from the Great08 RealNoiseKnown 
branch, for which the galaxies are modeled as either of bulge or 



0.04 ■ 
0.03 
0.02 
0.01 
0.00 ■ 

-0.01 - 
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0.2 
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-0.2 
-0.4 
-0.6 



1 2 3 4 5 5 

maximum order 

Figure 8. Offset a (top) and deviation from unity slope b — I 
(bottom) of GreatOS data sets given in Table [2] for different 
maximum orders, A^max, of sersiclet decomposition. Horizontal 
dashed lines indicate the ideal case of a = and 6=1. Errors of 
a and b were estimated by bootstrap fitting of f(x) = a + bx to 
the estimated and input shear values. 



However, the real lesson from Fig. [S] is that for A^max > 
the slopes are significantly below 1, i.e., the shear is signif- 
icantly underestimated. The higher-order sersiclet modes, 
which should account for any deviation the actual galaxy 
shape exhibits with respect to the axisymmetric Sersic pro- 
file, do more harm than good as they become increasingly 
prone to undersampling. Low resolution is a generic fea- 
ture of weak-lensing data, i.e., shear estimates based on 
sersiclet decompositions always suffer from substantial un- 
dersampling effects. Oversampling of the model within each 
pixel could cure this, but is - at least in our implementation 
- computationally infeasible. Finally, we note that we ob- 
served a bias in the shear estimates based on Sersiclets. The 
first component of the shear is systematically overestimated 
whereas the second component is systematically underesti- 
mated. We do not yet understand the precise origin of this 
bias but we may speculate that this is a pixellation effect 
(cf. Sect. 5.3 of |Massey et~ar||2007| ) . 



5 OUTLOOK: ORTHONORMALISING 

HIGHER-ORDER TAYLOR EXPANSIONS 

As discussed in Appendix [Xj the Sersic profile is the first- 
order Taylor expansion of any light profile. This naturally 
leads us to the expectation that with improving imaging 
quality the Sersic profile will not be a good match anymore, 
because it is "only" a first-order expansion. Consequently, 
an obvious strategy to enhance the Sersic profile is to allow 
for higher orders in the Taylor expansion of Eq. ( |A2[ ). Such 
higher-order radial profiles then can also be orthonormalised 
in order to describe azimuthal structures. This approach ap- 
pears to be very promising for investigations of galaxy mor- 
phologies, e.g., in the context of classification. However, it 



disc type, i.e. having Sersic index of ng 
tively. 



1 or 715 = 4, respec- 
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Figure 9. Comparing a Sersic profile witli ns = 2 (blue dasiied 
line) to a third-order profile with B = l/ng = 0.5, C = 0.1 and 
D = 0.25 (red sohd hnc). 



introduces further nonlinear model parameters which ren- 
ders it inappropriate for investigations of weak lensing where 
notoriously little data is available to constrain the model. 



5.1 Third-order profiles 



As explained in Appendix |A2| the next useful higher-order 
expansion which exhibits the correct boundary behaviour 
beyond the Sersic profile is a third-order expansion, 



p{r) ~ exp 



^A+Blog{r/l3)+Clog^{r/l3)+Dlog^(r/l3) 



] , (6) 



where D > and A, B, and C are arbitrary. Figure |9] shows 
an example of such a third-order profile in comparison to a 
normal Sersic profile. Evidently, the third-order profile can 
overcome two essential problems of pure Sersic profiles, since 
it (a) can exhibit a central cusp, which is also observed in 
real galaxies, and, (b) approaches zero for increasing radii 
faster than the pure Sersic profile. Therefore, we may specu- 
late that such higher-order Taylor expansions could provide 
a reasonable generalisation, if deviations from the normal 
Sersic profile are observed while azimuthal structures are 
still absent. For instances, this may help to describe the 
light profiles of elliptical galaxies or unbarred SO galaxies. 

In Fig. [To] we compare the performances of the third- 
order profile and the Sersic profile by fitting an elliptical 
galaxy from the SDSS database. Clearly, the residual map 
of the third-order profile is almost perfectly random noise 
whereas the residual map of the Sersic profile reveals system- 
atic mismodelling. Correspondingly, the ratio of x^-values 
of the third-order profile over Sersic fit is xi/Xs ~ 0.48. 
In fact, the third-order profile fits the data so well that we 
should ensure that it is not an overfit. For this purpose. 
Fig. [Tl] displays the distributions of normalised residuals for 
both models in comparison to the unit Gaussian (e.g. see 
[Andrae et al.|2010c| . Evidently, the normalised residuals of 
the Sersic profile have a broader distribution than the unit 
Gaussian which is indicative of underfitting the data. The 
normalised residuals of the third-order profile are closer to 
the unit Gaussian, i.e., they are closer to the truth. How- 
ever, their distribution does not peak sharper than the unit 
Gaussian, i.e., the third-order profile does not overfit the 
data. 




Figure 11. Model comparison via distributions of normalised 
residuals. The blue histogram shows the distribution of nor- 
malised residuals from the Sersic fit. The red histogram are 
the normalised residuals from the third-order profile. The black 
dashed line is a unit Gaussian. 



5.2 Numerical orthonormalisation 

Third-order profiles such as Eq. (|6| can be orthonormalised, 
too. Unfortunately, it is not possible to do this orthonormal- 
isation analytically]^ Therefore, the orthonormalisation has 
to be performed numerically. We start with the radial mono- 
mials (r^ ,r'' , . . .) which are linearly independent but 
not orthonormal. Using the definition of the scalar product 



according to Eq. (Bl I and Dirac notation, we first normalise 
the ground state. 



10) = 



(7) 



Second, we compute the first-order state 

|1) = (f-(0|f|0))|0) , 

which is then also normalised such that (1|1) — 1. Here, we 
used the following definition. 



(8) 



{k\r\l) = / drrBk{r)rBi{r) 



(9) 



where Bi denotes the radial basis function of order I. The 
basis functions of order I ^ 2 are then computed using the 
three-term recurrence relation 

10 = {f-{l- l\f\l - 1)) {I- l\f\l - 2)1/ - 2) . 

(10) 

This requires less numerical integrations than the brute- 
force Gram-Schmidt algorithm. Hence, it is faster and nu- 
merically more stable than the Gram-Schmidt algorithm. 

5.3 Revisiting the problems of sersiclets 

The conceptual problem of sersiclets was the postulated 
relation of steepness of the weight function and scale of 
polynomial oscillations that is not obeyed by real galaxies. 
This fixed relation is loosened by the third-order profiles. Of 
course, the Sersic index - or rather B — — - still has an 
influence on the oscillation scale of the radial polynomials. 
However, now there are two further model parameters C and 

^ The reason is that the substitution in Appcndix^is not bijec- 
tive anymore. 
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Figure 10. Fitting Sersic and third-order profile to an elliptical galaxy. The panels denote residual map of Sersic fit (a), Sersic model 
(b), original data (c), best-fitting third-order profile (d), and residual map from third-order profile (e). Panels, b, c, and d have identical 
scaling. The residual maps (panels a and e) both use plot ranges from —ba to 5a-. 




Figure 12. Undersampling test for third-order profiles. We com- 
pare basis sets of the Sersic (dashed lines) and third-order profile 
(solid lines) shown in Fig.|9] The test is performed on a 100x100 
pixel grid and the weight function has no intrinsic ellipticity. In 
contrast to Fig. [t] we have chosen 6 = 4— 1/3 in this case. Panel 
a: Determinant of coefHcient covariance matrix of ■ X which 
should be 1. Panel b: Largest and smallest diagonal elements of 
of • X which should be 1. Panel c: Largest absolute value of 
off-diagonal elements of of • X which should be 0. 



D which also influence the polynomials. Consequently, basis 
functions based on third-order profiles are more flexible and 
do not impose such a rigid relation as sersiclets do. 

The practical problem of sersiclets was undersampling, 
in particular for Sersic indices ns ^ 1. In fact, Fig. [9] al- 
ready suggests that third-order profiles can overcome this 
problem because they can exhibit central cusps instead of 
peaks. Therefore, the weights are not as highly concentrated 
and the polynomials may not oscillate as rapidly. We demon- 
strate this by repeating the orthonormality test of Fig. [T] for 
the third-order profile with parameters given in Fig. [orand 
its set of basis functions. Figure [12] shows the test results. 
In direct comparison to sersiclets with ns = = 2, the or- 
thonormalised third-order profile indeed suffers less strongly 
from undersampling, which is most obvious for the determi- 



nant of X"^ ■ X. In particular, there is a regime 5 5S /3 5C 15 
where undersampling is overcome and boundary truncation 
has not yet set in. Such a regime does not exist for the cor- 
responding set of sersiclets shown here. Consequently, third- 
order profile can indeed overcome the undersampling prob- 
lem of sersiclets. Nevertheless, there are parameter choices 
for third-order profiles where the undersampling problem is 
as bad as for sersiclets, e.g., when C = D = Q and the 
third-order profile reduces to a Sersic profile. However, by 
choosing appropriate priors it might be possible to exclude 
such parameter values in order to avoid undersampling. 



5.4 Computational feasibility 

We have shown that orthonormalisations of third-order pro- 
files can overcome the limitations of sersiclets while pre- 
serving their mathematical justification. However, we have 
not yet mentioned a serious practical limitation of this 
novel approach: The numerical orthonormalisation of third- 
order profiles is computationally highly expensive. Fitting a 
galaxy on a, say, 100x100 pixel grid while freely adjusting 
all model parameters is infeasible on a standard computer. A 
problematic work-around could be to first fit a simple third- 
order profile to the galaxy and then only orthonormalise 
the best-fit profile in order to model the data's deviation 
from the mean profile. However, this step-wise approach by 
construction is very unlikely to find the best-fitting basis- 
function expansion because the first step produces a biased 
fit and not a mean profile. Therefore, we do not recom- 
mend this approach. Another work-around would be to set 
up a library of basis functions for all realistic parameter 
values of third-order profiles, such that numerical orthonor- 
malisation has not to be conducted on-the fly during the 
fit. This approach is hampered by the fact that third-order 
profiles have three free parameters such that a decently de- 
tailed sampling of this three-dimensional parameter space - 
which is necessary in order to provide reliable fits - would 
produce an extensive library. The best solution is certainly 
to maintain the numerical orthonormalisation during the fit 
and to employ Graphics Processing Units (CPUs) instead of 
normal CPUs. For such a pure "number crunching" like nu- 
merical orthonormalisation, CPUs impressively outperform 
CPUs (e.g. [Fluke et al.|2011| ). 
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6 CONCLUSIONS 

In this article, we reinvestigated the method of sersiclets in 
order to overcome the limitations of shapelets reported by 
[Melc hior et al.| ( |2010[ ) in parametrising galaxy morphologies. 
Orthogonalising the Sersic profile is justified by the fact that 
the Sersic profile is the first-order Taylor expansion of any 
real light profile (Appendix |A1[ ). We obtain the following 
results: 

• From general radial light profiles like the Sersic profile, 
we can only construct polar basis function because Carte- 
sian basis functions do not exhibit the azimuthal symmetry 
required in the context of galaxy morphologies. Shapelets 
are the only exception to this rule. 

• Sersiclet basis functions exist analytically, such that no 
numerical orthonormalisation is required and models are 
simple to evaluate. Polar shapelets are a special case of the 
sersiclet basis functions. 

• Sersiclets outperform shapelets in the case of edge-on 
discs and elliptical galaxies, providing models with sub- 
stantially lower residuals. In particular, sersiclets overcome 
the ring-like artefacts introduced by shapelets when fitting 
galaxies with steep light profiles. However, in the case of ex- 
tended objects with lots of substructure, e.g., face-on spiral 
galaxies, shapelets are competitive. 

• Sersiclets indeed solve the problems of shapelets. How- 
ever, we revealed two new problems of this approach: 

- Sersiclets are prone to undersampling effects, if galax- 
ies are not well-resolved. Table [T] provides an estimate of 
the required resolution in order to avoid undersampling. 
This undersampling problem stems from the polar ba- 
sis functions being evaluated on a Cartesian pixel grid. 
The problem is therefore most prominent in the central 
regions of the model, where unfortunately also most of 
the galactic light is concentrated. Undersampling gives 
rise to an increased uncertainty of the model parameters. 
For decreasing resolution, increasing parameter uncertain- 
ties render sersiclet models less informative, which beyond 
some point obstructs most if not any application. If com- 
putationally feasible, the undersampling problem may be 
alleviated by oversampling the model, although oversam- 
pling inevitably compromises the interpretation of the ex- 
pansion coefficients. 

- Sersiclets postulate a relation between steepness of 
the weight function and spatial scale of polynomial oscilla- 



tions (cf. Sect. 2.3 1. However, real galaxy morphologies do 



not necessarily obey this relation, e.g., in the steepness of 
their bulge vs. the size and distribution of star-formation 
knots[^ Therefore, we have to expect modelling problems 
for sersiclets, which may not have a generally predictable 
impact on a given application. 

• We tested the performance of sersiclets by analysing 
simulated weak-lensing data from the GreatOS challenge 



(Bridle et al. 20091. We observed that higher-order modes 



do more harm than good. Given the typically low resolu- 
tion of realistic weak-lensing images and the undersampling 
problem of sersiclets, we have to conclude that sersiclets are 
inappropriate for shear measurements in weak lensing. 



In fact, this problem also applies to shapelets. 



We conclude that although the now formally correct sersiclet 
approach appears very reasonable at first glance, it suffers 
from substantial problems in practice. Therefore, for every 
usage of sersiclets, we recommend a critical assessment of 
the impact of undersampling and the postulated unphysical 
relation between steepness and oscillation scales onto the 
results. 

Finally, we demonstrated that third-order Taylor ex- 
pansions of galaxy light profiles can drastically improve the 
modelling fidelity in comparison to Sersic profiles. If or- 
thonormlised, third-order profiles can overcome the under- 
sampling problem of sersiclets and loosen the tight relation 
between steepness of the weight function and spatial scale 
of polynomial oscillations. Furthermore, they preserve the 
mathematical justification of sersiclets. However, such basis 
functions are computationally highly demanding. In fact, 
they may only be feasible using appropriate computer hard- 
ware such as GPUs. 



Jimenez-Teja & Bemtez (2011 1 introduce a set of basis 



functions based on Chebyshev rational functions. The au- 
thors demonstrate that their method is capable of providing 
excellent model reconstructions of galaxies of very differ- 
ent morphological types with very little computational cost. 
Orthonormalised higher-order profiles may provide similarly 
good models but certainly not at a similarly low computa- 
tional cost. Hence, the method of 'Jimenez-Tej a fc Bem'tez| 
(20111 appears to be superior]^ In general, basis-function 
expansions can describe azimuthal structures of galaxy mor- 
phologies very accurately. However, the astrophysical inter- 
pretation of the expansion coefficients is not obvious and has 
to be laboriously deduced. This is a generic problem of basis- 
function expansions. Conversely, Galfit3 ( Peng et al.|2010 l 
allows to directly model the azimuthal structures exhibited 
by galaxy morphologies, such as spiral-arm patterns. The 
modelling and choice of components are "guided" by astro- 
physical intuition. Consequently, results obtained from Gal- 
fit3 are easier to interpret. As discussed by | Jimenez- Tej a"^ 



Bem'tez (20111, this requires substantial manual interaction 



of the user which is conceptually hard to automatise. The 
above mentioned methods highlight a fundamental conflict 
for any attempt of modelling galaxies. One can either seek 
to first optimally describe their shapes and then interpret 
the abstract model parameters later - or attempt to imme- 
diately explain galaxy shapes as mixtures of a limited set 
of physically motived features. The sersiclets approach pro- 
vides a compromise between these two extremes, as it incor- 
porates the fiexibility and completeness of basis- function ap- 
proaches and the intuition of using a weight function, which 
is well-adapted to galactic morphologies. 

The sersiclet code is implemented in C-|--|- and is avail- 
able on request. 



Nevertheless, Chebyshev rational functions still have a long way 
to go. First, the choice of scale radius suggested by | Jimenez-Teja] 



I&: Benitez] | |2011| is a work-around that may erase discrimina- 
tive information when it comes to classifying galaxy morpholo- 
gies. Second, the authors yet have to investigate how their basis 
functions fare in the regime of poorly resolved galaxies which con- 
stitute the vast majority of galaxies in surveys. 
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APPENDIX A: SERSIC PROFILE AS 
FIRST-ORDER TAYLOR EXPANSION 

We reveal the nature of the Sersic profile as a first-order 
Taylor expansion. Our emphasis is on the justification of the 
choice of the Sersic profile for orthogonalisation]^ Further- 
more, we briefiy outline an alternative approach to generalise 
Sersic profiles by higher-order expansions. 

Al First-order expansion 

Let us consider the real, two-dimensional light profile of a 
galaxy I{x) projected onto the sky, which may exhibit arbi- 
trary radial and azimuthal structures. Due to observational 
effects, e.g., noise and resolution, the observed light profile 
Iobs{x) will be different. If noise and resolution have a strong 
impact, we cannot identify azimuthal structures anymore 
and only the radial decline is left, i.e., Iobs{x) « Iobs{r). Let 
us further consider a rescaling of the observed light profile 
p(r) — Iobs{r) / lobaiO), such that < p{r) ^ 1. Then we can 
take the logarithm of p{r) and due to logp(r) ^ 0, for r > 
we can also take the logarithm a second time, introducing 

p{r) = log(- logp(r)) . (Al) 

We now Taylor expand this function p{r) in log r at a char- 
acteristic radius log /J to first order, i.e., a constant plus the 
first nontrivial term, which yields, 

p(r) + B(logr - log/3) = A + log(r//3)^ . (A2) 

Let us transform backwards now, 

logp(r) = -e''<"> ^ -e'*(r//3)^ (A3) 

and hence 

p(r)«exp[-e'*(r//3)^] . (A4) 



^ [Ciotti fc Bertin] | |1999| also investigate a Taylor-expansion of 
the Sersic profile, but they expand in powers of Sersic index and 
not in powers of radius. 
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All we need to do now is to identify the coefficients A and B 
of the Taylor expansion in Eq. (A2|. If we simply "rename" 
these constants by i? = 1/ns and A = log(&„), then we have 
arrived at the definition of the Sersic profile 



ps(r) = exp 



-bn 



l/ns 



(A5) 



Evidently, the Sersic profile is a first-order Taylor expansion 
at ro = /3 > 0. In other words, in the limit of low resolu- 
tion and low signal-to-noise ratios (where substructures are 
negligible) any radial light profile is approximately a Sersic 
profilej^ 

This naturally explains why the Sersic profile is 
such a good fit to real galaxies in this regime. Furthermore, 
this also leads us to the expectation that with improving 
imaging quality the Sersic profile will not be a good match 
anymore, because it is "only" a first-order expansion. 



A2 Higher-order expansions 

Given the fact that the Sersic profile is the first-order Tay- 
lor expansion of a real light profile, an obvious strategy to 
enhance the Sersic profile is to allow for higher orders in the 
Taylor expansion of Eq. (A2l. However, a realistic profile 



has to be unity at r = and it has to approach zero for 
r — > oo. As the Taylor expansion is in logr, only expansions 
where the leading order term is of odd power in logr and 
has a positive expansion coefficient can satisfy these con- 
straints. Therefore, the next useful higher-order expansion 
beyond the Sersic profile is a third-order expansion, 

A + B log(r//3) -I- C log" (r//?) + D log^ (r//?) 



p{r) 



(A6) 



where D > and A, B, and C are arbitrary. 



APPENDIX B: MATHEMATICAL DERIVATION 
OF RADIAL PARTS OF SERSICLETS 

In this appendix, we give the derivation of the radial part 
of the sersiclets, showing that the basis functions have an 
analytic form. Our starting point is the scalar product of 
two radial modes of orders / and I' , 



no 

m = / 

Jo 



drrRi{r/P)Rv{r/l3) exp 



l/ns 



(Bl) 

where we have adopted Dirac notation from quantum me- 
chanics. Ri (r) denotes the radial polynomials we are looking 
for and the Sersic profile acts as the weight function or "met- 
ric" of this scalar product. Note the functional determinant 
drr, which is due to our integration in polar coordinates. 
We now change variables according to 



i{r) = 6 



1/t7 



such that 



, a2 ns 2nq 

drr = p -- — u 



^du . 



(B2) 



(B3) 



This is not true for radial profiles that are not difTerentiable, 
i.e., where no Taylor expansion exists. However, such profiles are 
unphysical. 



The limits of integration do not change under this transfor- 



mation. Then Eq. (Bl I reads 



(B4) 



The "new" weight function of this transformed scalar prod- 
uct is now of the form u*e~", where k = 2ns — 1, and the 
corresponding set of orthogonal polynomials are the associ- 
ated Laguerre polynomials^ 



d' 



U du' 



( -u l+k\ 



(B5) 



The L; exist if and only if fc > —1. This is guaranteed, since 
k = 2ns — 1 and the Sersic index satisfies ns > 0. The 
normalisation factor is 

' du Lf (u)Ll {u)u'^e-= + , (B6) 

where F denotes the Gamma function. Consequently, the 
radial parts of the sersiclets read 



Ri(r) 



oxp 



l/ns 



with normalisation factor 



Ni 



^ns r{l + 2ns) 



62" 



n 



(B7) 



(B8) 



Mind the factor of 1/2 in the exponent of Eq. (B9l, which 



arises from our definition that the Sersic profile is the weight 
function in Eq. ( |B1[ |. We could also have defined Eq. ( |B9[ ) 
with a pure Sersic profile, resulting in a factor 2 arising in 



Eq. (Bl I. Fit results obtained from both definitions will not 



differ, since these factors are absorbed in the parameter b. 

In the special case of polar shapelets these expressions 
simplify. First, for ns ~ 0.5, we note that k — 2ns — 1 = 0, 
which implies r(/-|-fc-|-l) = r(Z-f 1) = l\ and simplifies the as- 
sociated Laguerre polynomial Lf to the "normal" Laguerre 
polynomial = Li. Second, in order to obtain shapelets, 
we need to set fe = 1. Therefore, the radial parts of polar 
shapelets read 



1 



Riir) = Jg 



Xrlpy] exp 



2/32 



(B9) 



The basis functions introduced by 'Massoy fc Refregier| 
( ]2005| ) as "polar shapelets" (their Eq. (8)) differ from Eq. 
|B9|) by using a;'"''l/'^|^|j^2(x^) instead of Li{x^). Their 
basis functions are not orthogonal, e.g., consider 



{l = l,m = l\l = 2,m = 0) 



■ -r^ |1| r 1 /- 2\ |0| rO/ 2s 

drre r' L^lr )r' Li(r ). 



(BIO) 

which equals —^/tt/8 and is thus non-zero. Consequently, 
the "polar shapelets" defined by |Massey fc Refregier| ( [2005[ ) 
are incorrect. However, the linear independence is preserved 
because no two states with identical values of m (otherwise 
the azimuthal parts e""''' would differ) and different values 

2\ 



of I can have identical ji'^^'l'^i 



http:/ /mathworld. wolfram.com/LaguerrePolynomial. html 
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APPENDIX C: OPTIMISATION PROCEDURE 

Having defined sersiclet models, we also need to clarify how 
to fit such models to given imaging data. 



CI Maximum-likelihood solution 

For the moment, we assume that we are given estimates 
of the nonlinear model parameters {Nmax,ns, l3,xo,e). We 
define a vector c of free parameters that form the coefficients. 
Furthermore, we write the image as a vector y of size A'^, 
where A*' denotes the number of pixels in the given image. 
Finally, we define the so-called N x P design matrix X, 
such that Xnp is the p-th basis function evaluated at the 
n-th image pixel (after forward convolution with the PSF). 
With these definitions, we define the residual vector 



R^y-X- 



(CI) 



which is an vector of size A*' itself. Here "•" denotes matrix 
multiplication. As the pixel noise of modern imaging data is 
usually Gaussian in excellent approximation, we are allowed 
to define 



(C2) 



where E denotes the N x N pixel covariance matrix. 
In the case of uncorrelated Gaussian noise with constant 
variance in all pixels, E takes the simple form E — 
diag(a^, . . . , u^) — a^I. It is straightforward to show that 
the maximum-likelihood estimate that minimises is 

S ^ {X^ ■ T.'^ ■ X)~^ ■ X'^ ■ T.'^ ■ y . (C3) 

The matrix X"^ • E~^ ■ X has a special meaning: It is the 
P X P covariance matrix of the coefficients c, which would 
follow from a Fisher analysis. In this formalism, we get this 
information for free. 

The advantage of expressing the maximum-likelihood 
solution in terms of matrix operations is that we can employ 
fast and efficient linear-algebra algorithms. In our case, the 
basis functions are implemented in C++ and we use wrapper 
classes to import the packages ATLA^ and LAFACKp] 



Within these reasonable parameter ranges, we employ fiat 
priors. 

C3 Estimating the maximum order 

In the case of shapelets, the maximum order, A^max, is esti- 
mated via a reduced x^- However, demanding that equals 
the number of degrees of freedom is justified if and only if 



the model is purely linear (e.g., Barlow|[l993 Andrae et al 



2010c I. In the case of sersiclets, this assumption is definitely 
violated, since sersiclets contain many nonlinear fit parame- 
ters. We rather recommend to estimate A'^max by comparing 
the normalised residuals for models of different maximum or- 
ders. The optimal A'max is then defined by the model whose 
normalised residuals are closest to a Gaussian with mean 
zero and variance one ( Andrae et al.|2010c l. 

Concerning large samples of galajcy images, one can also 
decompose all objects using identical A'^max. This approach 
may be favourable because many techniques to analyse the 
resulting catalogue of models require that all coefficient vec- 
tors have the same dimensionality (e.g. clustering analysis, 
cf. [Kelly fc McKay||2004) ,2005) |Andrae et al.|2010b l. 



This paper has been typeset from a TjnjX/ IWJJC file prepared 
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C2 Optimising the nonlinear parameters for given 

A^ max 

We now assume that we are given an esimate of the maxi- 
mum order, A'max, and need to get estimates of (ng, j3, xo, e). 
These estimates are derived from a Simplex algorithm 
( Nelder fc Mead||1965l ) that we incorporate from the GNU 
Scientific Library (GSLp] The Simplex algorithm is an it- 
erative optimisation algorithm that does not employ deriva- 
tives but only evaluations of x^- We also employ priors in 
order to restrict the model parameters to reasonable values. 
These priors in detail ensure that: 

• The Sersic index is in the range 0.1 ^ ns ^ 8. 

• The scale radius satisfies /3 > 0. 

• The object centroid xq is within the pixel grid. 

• The complex ellipticity satisfies jej < 1. 



http: / / math-atlas.sourceforge.net 

http: / /www. netlib.org/lapack/ 

http: / /www. gnu.org/software/gsl/manual/ 



